home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / endo / north.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-03  |  5.9 KB  |  220 lines
  1. /*************************************************************************
  2.  *                                                                       *
  3.  *  Copyright (c) 1992, 1993 Ronald Joe Record                           *
  4.  *                                                                       *
  5.  *  All rights reserved. No part of this program or publication may be   *
  6.  *  reproduced, transmitted, transcribed, stored in a retrieval system,  *
  7.  *  or translated into any language or computer language, in any form or *
  8.  *  by any means, electronic, mechanical, magnetic, optical, chemical,   *
  9.  *  biological, or otherwise, without the prior written permission of:   *
  10.  *                                                                       *
  11.  *      Ronald Joe Record (408) 458-3718                                 *
  12.  *      212 Owen St., Santa Cruz, California 95062 USA                   *
  13.  *                                                                       *
  14.  *************************************************************************/
  15.  
  16. #include <math.h>
  17. #ifndef NeXT
  18. #include <values.h>
  19. #endif
  20. #include "defines.h"
  21.  
  22. pair GNP, RK;
  23. double pB, A2T, CT, A1T, ET, FT, GT; 
  24.  
  25. #ifdef NorthSouth
  26.  
  27. pair
  28. north(x, y, p) /* North-South model from Chichilnisky, Abraham & Record */
  29. double x, y; 
  30. double *p;
  31. {
  32.     static pair coord;
  33.     extern void gamma_quadratic();
  34.     extern void compute_gnp();
  35.  
  36.     gamma_quadratic(x, y, p);
  37.     compute_gnp(x, y, p);
  38. #ifdef RESTRICT
  39. /*    Restrict to non-negative GNP */
  40.      coord.x = Max(0.0,(p[13] * GNP.x) + ((1.0 - p[12]) * x));
  41.      coord.y = Max(0.0,(p[1] * GNP.y) + ((1.0 - p[0]) * y));
  42. #else
  43. /* this does not restrict to non-negative gnp  */
  44.     coord.x = (p[13] * GNP.x) + ((1.0 - p[12]) * x);
  45.      coord.y = (p[1] * GNP.y) + ((1.0 - p[0]) * y);
  46. #endif
  47.     return(coord);
  48. }
  49.  
  50. pair
  51. dnorth(x, y, p)
  52. double x, y;
  53. double *p;
  54. {
  55.     extern void gamma_quadratic(); 
  56.     extern pair dmodels();
  57.  
  58.     gamma_quadratic(x, y, p);
  59.     return(dmodels(x, y, p));
  60. }
  61.  
  62. pair
  63. dmodels(x, y, p)
  64. double x, y;
  65. double *p;
  66. {
  67.     static double r, dpBdx, dpBdy, dxdx, dxdy, dydx, dydy, dnsq, dssq;
  68.     static double a2n, dn, c1n, lbarn, alphan, c2n, sn, deltan, cnsq, cssq;
  69.     static double a2s, ds, c1s, lbars, alphas, c2s, ss, deltas, a1n, a1s;
  70.     static pair z;
  71.  
  72. printf("");
  73.     z.x = z.y = MAXDOUBLE;
  74.     a1n = p[14]; a2n=p[15]; dn=p[18]; c1n=p[16]; lbarn=p[20]; 
  75.     alphan=p[21]; c2n=p[17]; sn=p[13]; deltan=p[12];
  76.     a1n = p[2]; a2s=p[3]; ds=p[6]; c1s=p[4]; lbars=p[8]; 
  77.     alphas=p[9]; c2s=p[5]; ss=p[1]; deltas=p[0];
  78.     dnsq = dn * dn; dssq = ds * ds;
  79.     cnsq = c1n * c1n; cssq = c1s * c1s;
  80.     if ((pB == 0.0)||(dn == 0.0)||(ds == 0.0)||(dnsq == 0.0)||(dssq == 0.0))
  81.         return(z);
  82.     r =  (1.0/dn) * ((c1n * lbarn) + (alphan * c1n * c2n / dn));
  83.     r += (1.0/ds) * ((c1s * lbars) + (alphas * cssq / ds));
  84.     r /= (a2n * x / dn) + (a2s * y / ds);
  85.     dpBdx = r * a2n / (2.0 * dn);
  86.     dpBdy = r * a2s / (2.0 * ds);
  87.     r = (alphan*c2n*c2n/dnsq) + (c2n*lbarn/dn) - (a2n*x/dn);
  88.     dxdx=(-1.0*pB*a2n/dn)+(r*dpBdx)+(a1n/dn)-((alphan*cnsq)/(dnsq*pB*pB));
  89.     dxdx *= sn;
  90.     dxdx += 1.0 - deltan;
  91.     dxdy = sn * ((r - (alphan*cnsq/(dnsq*pB*pB))) * dpBdy);
  92.     r = (alphas*c2s*c2s/dssq) + (c2s*lbars/ds) - (a2s*y/ds);
  93.     dydy=(-1.0*pB*a2s/ds)+(r*dpBdy)+(a1s/ds)-((alphas*cssq)/(dssq*pB*pB));
  94.     dydy *= ss;
  95.     dydy += 1.0 - deltas;
  96.     dydx = ss * ((r - (alphas*cssq/(dssq*pB*pB))) * dpBdx);
  97.     z.x = (dxdx * dydy) - (dxdy * dydx);
  98.     z.y = dxdx + dydy;
  99.     return(z);
  100. }
  101.  
  102. void
  103. solve_quadratic(x, y, p)
  104. double x, y;
  105. double *p;
  106. {
  107.     static double pB1, pB2, a, b, c, q;
  108.     extern double sqrt();
  109.     extern double   VT;
  110.  
  111.     a = (x*p[15]/p[18]) + (y*p[3]/p[6]);
  112.     b = VT * (-1.0); 
  113.     c = (p[9]*p[4]*p[4]/(p[6]*p[6])) + (p[21]*p[16]*p[16]/(p[18]*p[18]));
  114.     if (b < 0)
  115.         q = -0.5 * (b - (sqrt((b*b) - (4.0*a*c))));
  116.     else
  117.         q = -0.5 * (b + (sqrt((b*b) - (4.0*a*c))));
  118.     if (a == 0.0)
  119.         pB1 = 0.0;
  120.     else
  121.         pB1 = q/a;
  122.     if (q == 0.0)
  123.         pB2 = 0.0;
  124.     else
  125.         pB2 = c/q;
  126.     if ((pB = Min(pB1, pB2)) <= 0)
  127.         pB = Max(pB1, pB2);
  128. }
  129.  
  130. void
  131. gamma_quadratic(x, y, p)
  132. double x, y;
  133. double *p;
  134. {
  135.     static double pB1, pB2, a, b, c, q, De_sq, n_De_sq;
  136.     extern double sqrt();
  137.     extern double   WT;
  138.  
  139.     n_De_sq = p[18] * p[18]; De_sq = p[6] * p[6];
  140.     a = (p[21]*p[17]*p[17]/(n_De_sq))-(((x*p[15])+(p[17]*p[20]))/p[18]);
  141.     b = a * (1.0 - p[23]);
  142.     c = (p[9]*p[5]*p[5]/(De_sq))-(((y*p[3])+(p[5]*p[8]))/p[6]);
  143.     q = c * (1.0 - p[11]);
  144.     ET = b + q;
  145.     CT = ((1.0/p[18])*((p[16]*p[20])+(p[21]*p[16]*p[17]/p[18])-(p[14]*x)))+
  146.          ((1.0/p[6])*((p[4]*p[8]) + (p[9]*p[4]*p[5]/p[6]) - (p[2]*y)));
  147.     a = 2.0*p[21]*p[16]*p[17]/n_De_sq;
  148.     b = x*p[14]/p[18];
  149.     c = p[16]*p[20]/p[18];
  150.     q = 1.0 - p[23];
  151.     GT = q * (b - a - c);
  152.     a = 2.0*p[9]*p[4]*p[5]/De_sq;
  153.     b = y*p[2]/p[6];
  154.     c = p[4]*p[8]/p[6];
  155.     q = 1.0 - p[11];
  156.     FT = GT + (q * (b - a - c));
  157.     a = (1.0 - p[23]) * (p[21]*p[16]*p[16]/n_De_sq);
  158.     b = (1.0 - p[11]) * (p[9]*p[4]*p[4]/De_sq);
  159.     GT = a + b;
  160.     a = ET; b = CT + FT; c = GT - WT;
  161.     if (b < 0)
  162.         q = -0.5 * (b - (sqrt((b*b) - (4.0*a*c))));
  163.     else
  164.         q = -0.5 * (b + (sqrt((b*b) - (4.0*a*c))));
  165.     if (a == 0.0)
  166.         pB1 = 0.0;
  167.     else
  168.         pB1 = q/a;
  169.     if (q == 0.0)
  170.         pB2 = 0.0;
  171.     else
  172.         pB2 = c/q;
  173. /*    if ((pB = Min(pB1, pB2)) <= 0)  */
  174.         pB = Max(pB1, pB2);
  175. }
  176.  
  177. void 
  178. compute_rk(x, y, p)
  179. double x, y;
  180. double *p;
  181. {
  182.     if (p[18] != 0.0)
  183.         RK.x = (p[14] - (p[15]*pB))/p[18];
  184.     else {
  185.         printf("Setting North RK to MAXFLOAT\n");
  186.         RK.x = MAXFLOAT;
  187.     }
  188.     if (p[6] != 0.0)
  189.         RK.x = (p[2] - (p[3]*pB))/p[6];
  190.     else {
  191.         printf("Setting South RK to MAXFLOAT\n");
  192.         RK.y = MAXFLOAT;
  193.     }
  194. }
  195.  
  196. void 
  197. compute_gnp(x, y, p)
  198. double x, y;
  199. double *p;
  200. {
  201.     static double n_L, s_L, Dsq;
  202.  
  203.     if (pB != 0.0) {
  204.         Dsq = p[18] * p[18];
  205.         n_L = (p[21]*p[17]*p[17]/Dsq) + (p[17]*p[20]/p[18]) - (p[15]*x/p[18]);
  206.         s_L = (-2.0*p[21]*p[16]*p[17]/Dsq)+(p[14]*x/p[18])-(p[16]*p[20]/p[18]);
  207.         GNP.x = (pB*n_L) + s_L + ((p[21]*p[16]*p[16])/(Dsq*pB));
  208.         Dsq = p[6] * p[6];
  209.         n_L = (p[9]*p[5]*p[5]/Dsq) + (p[5]*p[8]/p[6]) - (p[3]*y/p[6]);
  210.         s_L = (-2.0*p[9]*p[4]*p[5]/Dsq)+(p[2]*y/p[6])-(p[4]*p[8]/p[6]);
  211.         GNP.y = (pB*n_L) + s_L + ((p[9]*p[4]*p[4])/(Dsq*pB));
  212.     }
  213.     else {
  214.         printf("Setting GNP's to MAXFLOAT\n");
  215.         GNP.x = MAXFLOAT;
  216.         GNP.y = MAXFLOAT;
  217.     }
  218. }
  219. #endif /* NorthSouth */
  220.